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Abstract 

This paper presents a reliable algorithm for the evaluation of a quadratic perfor- 
mance index and its gradients with respect to the controller design parameters^ The 
algorithm is part of a design algorithm for optimal linear dynamic output-feedback 
controller that minimizes a finite-time quadratic performance index. 1 lie numerical 
scheme is particularly robust when it is applied to the control-law synthesis for sys- 
tems with densely packed modes and where there is a high likelihood of encountering 
degeneracies in the closed-loop eigensystem. This approach through the use of an ac- 
curate Pade series approximation does not require the closed-loop system matrix to 
be diagonalizable. The algorithm has been included in a control design package for 
optimal robust low-order controllers. Usefulness of the proposed numerical algorithm 
has been demonstrated using numerous practical design cases where degeneracies occur 
frequently in the closed-loop system under an arbitrary controller design initialization 
and during the numerical search. 


1. Introduction 

Traditional design methods in linear optimal control for continuous-time systems have been 
extensively treated in recent literature [1], Development of these control systems are usually- 
based on characterization of the control problem under the setting of optimization of the 
two-norm of a set of controlled output responses to random disturbance inputs or initial con- 
ditions. Additional consideration of design robustness is taken by formulating the problem to 
include Jf~- norm bound constraints for a class of additive and multiplicative uncertainties 
applied at the plant inputs and/or outputs. Solutions are obtained for both the state- an 
output-feedback design problems and involve in the majority of cases solving an appropriate 
set of algebraic Riccati equations [2,3]. Theoretical studies of these approaches have been 

-The work of B. VanSteenwyk and U. Ly is supported in part by NASA Ames Research Center under 
grant contract NAG-2-691. 


1 


the major concern of researchers in the control field and major breakthrough has been made 
in recent work by Stoorvogel [4,5]. An alternate and seldomly mentioned design option for 
robust multivariable control of linear time- invariant systems is based on direct numerical 
optimization of a quadratic performance index with an arbitrarily specified controller struc- 
ture. We believe that careful formulation of the design problem under nonlinear constrained 
optimization can be of great value in the synthesis of robust multivariable control systems. 

Early work in this area have been published by Levine and Athans [6], Anderson and 
Moore [7], and extensive review of the subject was performed by Makila [8]. Recently, a 
new look into parameter optimization to multivariable control synthesis is provided by Ly 
[9] where he used a quadratic performance based on finite-time horizon. In the latter work, a 
numerical optimization technique based on [10] was used. At each design iteration it requires 
a precise evaluation of a finite-time quadratic performance index and its gradients with 
respect to the design parameters. Analytical expressions have been developed to evaluate 
these quantities under the key assumption that the closed-loop system being diagonalizable. 
This assumption is found to be unsatisfactory and is the cause of convergence difficulties in 
the iterative search when it attempts to invert an ill-conditionned eigenvector matrix for the 
diagonalization. The work presented in this paper is to resolve this numerical difficulty and 
thereby extends the results of Ly [9] for cases where the closed-loop systems are degenerate, 

1. e the closed-loop system has repeated eigenvalues and the corresponding set of system 
eigenvectors does not span the whole state-space of the closed-loop system. 

The paper is organized as follows. Section 2 reviews the problem formulation for a linear 
optimal control design using direct parameter optimization. Analytical expressions for the 
evaluation of the quadratic cost function and its gradients with respect to the controller 
design parameters are also given in section 2. The current approach to evaluate these quan- 
tities are briefly reviewed in section 3. An alternate numerical scheme for the exponential 
matrix and convolution integrals involving exponential matrices is presented in section 4. 
Approximation methods for the evaluation of these matrix quantities using Pade series are 
described in details in sections 5, 6, and 7. A design algorithm based on these numerical 
solution schemes has been incorporated into a computer-aided design package described in 
[9], A simple design problem to motivate the need for a numerical algorithm that handles 
degeneracies in the closed- loop system matrix is given in section 8. Optimal solutions are ob- 
tained using the proposed method and the diagonalization method from [9]. The numeiical 
algorithm has also been applied to the synthesis of an active control system foi the JPL laige 
space structure developed under the LSCL research program [16]. Results of this application 
are presented in section 9. Conclusions are given in section 10. 

2. Problem Formulation 

In this section, we recall the problem formulation described in Ly [9] for the control synthesis 
of a robust low-order controller in a linear time-invariant system. The system P'(s) is 
controlled by a constant-gain controller C(s) as depicted in Figure 1 where ,t/*(s) is the 
controlled output vector, y' s {s) the measurement output vector, u>'(s) the disturbance input 
vector and u‘(s) the control input vector. As a consistent notation, the superscript i is used 
throughout this paper to denote the system model at the i th plant condition. Note that the 
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controller C(s) is considered to be fixed , i.e does not vary with the design condition. It is 
modelled as a linear time- invariant system of arbitrary order whose formulation accomodates 
both a feedforward and a feedback controller structures. Robustness requirement in the 
context of our problem formulation is defined under the conditions that the control system 
C(s) stabilizes the plant P‘(s) over a class of design conditions ( i = l,N p ). 

State equations describing the system model P'{s) of Figure 1 are as follows. Notice 
that, in the problem formulation, we assume without loss of generalities that all the system 
states are initially acquiescent. This assumption is not restrictive since one can always 
establish impulsive inputs tu'(t) together with the appropriate influence matrix to represent 
any state initial conditions. At the i th plant condition, the system design model is described 
by equations (l)-(3) below. 

State Equations: 

fp(t) = F‘x t {t) + G t u^t) + rw l (t) (1) 

1 x'(0) = 0 

where x'{t) is a n X 1 plant state vector, u'(t) an rn x 1 control vector, w'{t) an m' X 1 
disturbance-input vector, F’ an n x n state matrix, G l an n X m control distribution matrix 
and P an n x m! input-disturbance distribution matrix. 

Measurement Equations: 

y' s (t) = Hlx l (t) + D' $u ii'(t) + D' sw w l (t) (2) 

where y l s (t) is a p x 1 measurement vector, H\ a p x n state- to-measurement distribution 
matrix, D l su a p x m control-to-measurement distribution matrix and D sw a p x m input- 
disturbance- to- measurement distribution matrix. 

Criterion Equations: 

% /’(*) = Hlxft) + D l cu u\t) + D‘ cu ,w\t) (3) 

where y‘(<) is a p 1 x 1 criterion vector, H' c a p 1 x n state-to-criterion distribution matiix, 
Di u a p' x m control-to-criterion distribution matrix and D‘ cw a p' x in' input-disturbancc- 
to-criterion distribution matrix. 

For generality, the disturbances w*(t) are modeled as outputs ol a linear time- invariant 
system excited by either impulse inputs or white noises. In this manner, one can shape the 
disturbance signals to have any deterministic responses (e.g filtered step functions, sinusoidal 
functions, exponentially decayed or growing sinusoidal functions, etc...) or to model stochas- 
tic inputs with any given power spectral density functions. At the i th plant condition, the 
disturbance model is given by equations ( 4 ) - ( 5 ) below. 

Disturbance State Equations: 

fiL(0 = K<(i) + rw(t) (4) 

1 4(0) = o 

where x\ (/) is a n‘ x 1 disturbance state vector, rf{t) a m' x 1 vector of either parameterized 
random Impulses (i.e if(t) = t f 0 6(t) with E^} = 0, and E[ V 'X T ) = W 0 ), or white-noise 
processes i f(t) with zero mean and covariance )r/ lT (r)] = W 0 S(t — r). The matrix W a 

is an m' x m! diagonal positive semi-definite matrix, F' u , an ??' x n’ state matrix of the 
disturbance model and T^, an n 1 x in' input-distribution matrix. 
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Disturbance Output Equations: 

w‘(t) = Hyjt)+ (5) 

where w'(t) is a m' x 1 disturbance output vector, H' w an w! x n' disturbance output matrix 
and D l an m! x m' direct feedthrough distribution matrix. 

State model of the controller C(a) in Figure 1 is that of a linear time- invariant system 

described by equations (6)- (7) below. 

Controller State Equations: 

f z(t) = Az(t) + Byl(t) ( 6 ) 

U(0) = 0 

where z(t) is a r x 1 controller state vector, A a r x r state matrix of the controller and B 
a r x p measurement-input distribution matrix. 

Control Equations: 

u‘(t) = Cz(t) + Dy\(t) (7) 

where u'(t) is an m x 1 feedback control vector, C an m x r control-output distribution 

matrix and D an m x p direct feedthrough matrix. 

For control-law synthesis, all the elements of the controller state matrices can be chosen 
as design parameters and some of them can be left fixed at pre-assigned values. In addi- 
tion, if needed, linear and nonlinear equality or inequality constraints can be established 
among the selected design parameters in order to ensure a particular design structure. Foi 
convenience in the derivation of the performance index and its gradients with respect to 
the controller design parameters, we define a matrix C a that groups all the controller state 
matrices (A, 5, C, D) in one compact form as follows, 


C 0 = 


D C' 
B A . 


(m+r)x 


(p+ r ) 


( 8 ) 


Thus, specification of the matrix C 0 will completely define the controller state model. Ob- 
viously, for the case of a static output-feedback design (i.e the controller order r = 0), we 
simply have C Q = D. In Section 8, we will formulate a control design problem based on the 
minimization of a performance index using the controller defined in equations (6)-(7). 

To examine the entire class ol H ^- norm based control problems and to handle the piob- 
lem of sensitivity to plant modeling uncertainties, we define the objective function given in 
equations (9) and (10). This formulation turns out to be versatile and well-posed for the 
setting of a nonlinear constrained optimization problem. However, depending on the types 
of disturbance model, that is whether the disturbance outputs w‘(t) are responses to impulse 
or white-noise inputs, different definitions for the objective function are needed. They aie 
given below. 

(a) For random impulse inputs »/'(<): 


.7(1/ ) = 5 Z ( W i /o' e 2 "' 1 

l— 1 

E [ y' c T (t)Qy c {t) + u' T (t)R l u l {t) ] dt } 
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The expectation operator E [- ] is over the ensemble of the random variables if 0 in the param- 
eterized impulse inputs r/‘(f) = Control design problems formulated with the above 

performance index J(t f ) are often classified under the category of deterministic control. 
Under this category are, for example, the familiar control problems of command tracking 
control, disturbance rejection of unwanted but known external input signals, implicit and ex- 
plicit model-following designs, H 2 - control to initial conditions and //'^-control to sinusoidal 

inputs. 

(b) For random white-noise inputs 


I 






( 10 ) 


The expectation operator E a i[-] is over the ensemble of the random processes defined in the 
input variables t;‘( t ) for a closed-loop system destabilized by a factor a*. The destabilization 
effectively adds a value a* to the diagonal elements of the closed-loop system matrix. With 
the given performance index, one can address the entire class of H 2 - norm based control 
design problems. For example, we can solve for the linear quadratic regulator design (LQ), 
linear quadratic gaussian (LQG) design, loop transfer recovery (LTR), closed-loop transfer 
recovery (CLTR), model reduction based on a minimization of H 2 - norm of the error. 

Note that the performance indices given in equations (9) and (10) are evaluated to a 
finite- time horizon tj. The use of a finite time plays a significant role in the implementation 
of a reliable design algorithm for the optimum steady-state solution. It should be recognized 
that the objective function is well-defined regardless of whether the feedback control-law 
is stabilizing or not. Furthermore, a class of problems associated with command tracking 
of neutrally stable or unstable target responses (e.g step and ramp commands, sinusoidal 
trajectories) are only tractable under the setting of a finite-time objective function but not 
in the confine of a steady-state objective function where tj — > oo. In practice, steady-state 
results, whenever possible, are usually achieved when the terminal time tj is equal to five oi 
six times the slowest time constant in the closed-loop system responses. 

There are other unique features, besides the concept of design to a finite terminal time t f , 
that we have incorporated into the design objective function of equations (9)- (10). First ot 
all, this objective function is not the usual quadratic cost function defined in traditional linear 
optimal control problems. It is instead a weighted average of quadratic performance indices 
evaluated over the entire set of design conditions ( i = 1 ,N P ). Different weights are assigned 
to each plant condition through the scalar variable W* where W l p > 0. Of course, if N p = 1, 
then we recover the usual quadratic cost function for a single nominal design condition. 
The time-weighted factor e 20 ' 1 further allows us to impose directly a stability requirement 
for the closed-loop eigenvalues. Namely, when a steady-state design has been achieved and 
the optimum objective function is bounded, then closed-loop system eigenvalues for the 
controllable modes will have real parts less than -a 1 . Finally, the weighting matrices Q‘ and 
R l are symmetric and positive se? 7 ii- definite matrices. Note that our solution approach to the 
minimization of the objective function J(tj) based on nonlinear optimization does not require 
the control weighting matrix R l to be positive definite. In fact, in some design problems 
such as command tracking and model reduction, an objective function representing simply 
the tracking or model-matching errors does not include the control term, hence R‘ = 0. 


In this section, we provide analytical expressions for the evaluation of the objective 
function J(tj) and its gradients dJ/dC a with respect to the controller matrix C 0 . Details of 
the derivation can be found in [9]. For simple technicality, the problem formulation assumes 
that there is no implicit-loop paths within the feedback control system. Namely, the control 
input u’(t) or the measurement output y l s (t ) should not have any direct link to itself. This 
translates into the conditions that one of the product DD l su or D l au D must be zero. This is not 
restrictive since in practice either actuation or sensor dynamics would be incorpoiated into 
the design models and thereby resulting in a system that satisfies the above assumption, 01 
one can simply reformulate an equivalent problem with a set of measurement outputs where 
D' su = 0. Let’s assume without loss of generality that we have the case of DD' SU = 0, then 
the closed-loop system can be written in the form [9] shown in figure 2 or simply, 


i*(<) = F*x H {t) + r y (0 with a.- ,! (0) = 0 


where 


vil 

J (n-j-r+n') X 

(n+r+n') 

= 

F'+ 

G'C 

(TH 

G'DH\ 


G'DD' SW )H X W 

B{ 1 + 

A+ 

B(l+ 

D' su D)H' s 

BC' SU C 

Di u D)D' sw H' w 

0 

0 

K 


n 


n+r+n')x m* 


(P + G‘DD' s JDl 
B(l+Dl u D)D' sw D‘ u 


[F s } p x(n+r+n’) ~ 

[(I + D' su D)Hl D' SU C (I + D' su D)D' sw Hi) 

\d’;u »■=[(! + d:,d)D‘ iw di} 

[H c ]p'x(n+ r+n’) 

[HI + D' CU DH' S D' cu C (D' cu DDl w + D' CW )H' W 


and 


[C"U ,„ +r+ »., - [DHi C DDlX] 

With these definitions, equations (2), (3) and (7) for y % s {t),y' c (t) and u l {t) become 


( 11 ) 


( 12 ) 


( 13 ) 

(14) 

(15) 

(16) 

(17) 


vi(t) = H^(t) + D'; v \t) (18) 

y*(0 = H?x H (t) + {D i cu DDsw i + DiJDijfit) (19) 

u i (t) = C , 'x , '{t) + DD\ w D' w r i '(t) (20) 

For a well-posed performance index product of direct feedthrough term (D x cu DD' sw + 

D z cw )D' w in the criterion output y' c (t) and the penalty weighting matrix Q l must be zero. 
Similarly, product of the direct feedthrough term DD l sw D‘ w in the control output u‘{t) and the 
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penalty weighting matrix R l must also be zero. Under these circumstances, the performance 
index J{tj) in equations (9) and (10) can be written as, 


1 




( 21 ] 




where 


— fb AF'+o'I)t 


(22) 


L'(t f ) =fo'e { 

[H?Q'H l c + C iT R i C i ]el F ' +a ' I)Tt dt 

In the derivation of the analytical gradients of the performance index J(t } ) with respect 
to design parameters in C a of the controller state matrices, it is convenient to express the 
closed-loop system matrices in terms oi C 0 explicitly, as suggested in [9], 


where 


P + (G' 0 + T 2 C 0 D\)C 0 H' 0 

(23) 

r i + (G' o + T 2 C 0 D[)C 0 D' 0 

(24) 

H’^ = H[ + D\C 0 H[ 

(25) 

= [DH\ C DD\ w Hl\ 

(26) 


- T x C 0 Hi 

[^o](n+r+n')x (n+r + n') ” 

(m+r) ' 

[rj(7i+r+n')x m ‘ ““ 
wx+ r)x (n+r+n') 


p o r ! //‘ 
ooo 
0 0 F' w 

G' 0 
0 I 
0 0 

r D i 

o 

r:.. 


Hi o DlJll 
0 I 0 


[H\} p , x {n+r+n)) = [Hi 0 D' CW H' W ] 

[^o](p+>9x m' = 


D l D l 

sw w 


[^l](p+»-)x (m+r) — 


D'su 0 
0 0 


[^]p'X (m+r) — [D cu 0] 


(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 
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( 35 ) 


[7l]mX (m+r) [I 0] 

0 0 
0 I 
0 0 


[72](n+r+n')X (m+r) 


P 3 ] (p+r') X (n+r+n') — 


0 0 0 
0 I 0 


(36) 


(37) 


It has been shown in [9] that derivative of the performance index J(t f ) with respect to 
the controller matrix C 0 (i.e dj fdC 0 ) can be obtained explicitly from the following set of 

equations, 

dJldCo = 

T&I + Tj T RC' i )X i {tj)H™+ (38) 

(G l 0 + T 2 C 0 D\) T [M l {t s )H^ + C{t f )T ,t W^}+ 

Tj[M'(i f )H? + C'(tj)r’'W' 0 D?}(D\C 0 ) T }} 


dt 


(39) 


(40) 


(41) 


where 

*•(*/) = 

jtf e (F''+al)t p/t|^/ip^T e (F /l +aI )' 1 
£*(</) = 

/'/ £ (F" ta[) T « ; //"’I Q //'' + C ,a h,l) ' dt 

MV,) = /o'/o e' F " + °W-°>(H?QHZ + 

C" T RC' l V r " + " l|< p Wl P 1 F " 1 ''''' " t/a '// 

From equations (21)-(22) and (38)-(41), evaluation of the performance index J(t } ) and 
its gradients dJ/dC 0 would require algorithms for efficient computation of the follow- 
ing two types of integrals of matrix exponentials: X(t) = f 0 e Ar Be Cr dr and M{t) — 

io Jo Be Cv De Es ds dv. In the next section, we review briefly the numerical algorithm 

developed in [9] for the computation of ■ V ( / ) and AA(t). 

3. Current Method for Evaluating X(t) and M{t) 

Previous methods for evaluating X{t) and M(t) involve basically the diagonalization of the 
matrices A, C and E in the exponential functions. The procedure requires the determination 
of the eigenvalues and eigenvectors of these matrices. It is further assumed for convenience 
that similarity transformations can be constructed from these eigenvectors to diagonalize 
the respective matrices. Namely, there exist nonsingular transformations V A and V c and V E 
such that 

A = V a AaYa 1 , C = V C A C V C 7 1 , E = V E \ E V E (42) 

where the matrices A 4 . A q and A /-; are diagonal. Under these assumptions one can expiess 
for example the exponential function of e At as 


■A t - 6 v a^aV a ‘t _ v A e AAt V y ~ l 


(43) 


8 


Usage of this decomposition in the calculation of X(t) is shown below. 

X(t) = fe AT Be CT dT 
Jo 

= V A {f\ AAT Be AcT drJ VJ 1 (44) 

where B = V^ 1 BVc- Advantage of this approach is based on the fact that the exponential 
function of a diagonal matrix is also diagonal. In this case, time integration in X{t) can be 
performed directly by explicit integration of product of scalar exponential functions. The 
resulting numerical algorithm is quite accurate and efficient, provided that the transfoima- 
tion matrices V A and V c are not ill-conditionned. A similar procedure can also be applied 
to the evaluation of M{t). Complete discussion can be found in Appendices C and D of 
[9], However, breakdown of this algorithm will occur when the matrices A , C or E become 
degenerate or near degenerate; a situation that becomes eminent when we address control 
of flexible structures with densely packed modes as demonstrated in the design examples of 
sections 8 and 9. 

Clearly, in order to have a reliable design algorithm for optimal low-order output-feedback 
control synthesis [9], one must develop a robust numerical scheme to evaluate matrix integrals 
of the form shown in X(t) and M{t) for the case of a degenerate system. 


4. Alternative Approaches for Solving X(i) and M(t) 

One rather simple approach is to evaluate 

X(t) = f e Ar Be CT dr , (45) 

Jo 

M(t)= f / U e A{v ~ 3) Be Cv De Es ds dv (46) 

Jo Jo 

directly using techniques based on numerical quadrature. Efficiency of numeiical integration 
techniques is poor; especially when it requires small integration step size foi satisfactory 
accuracy in the case of stiff system matrices A, C and E. Another possibility is to use some 
types of algebraic Lyapunov equations for the solution of X(t) and A4(f). For example, it 
can be easily shown that the matrix X{t) can be obtained from the solution of the following 
Lyapunov equation, 

AX{t) + X(t)C = [e AT Be CT ] Q (47) 

Solution of equation (47) exists if A,(A) + Aj(C) ^ 0. This condition will not be satisfied 
in general for arbitrary system matrices A and C . Thus, from practical purposes X(t) and 
likewise M{t) cannot be solved from a scheme based on Lyapunov equations. 

Another possible approach is based on the direct use of exponential matrix. It is well- 
known [12] that convolution integrals involving matrix exponentials, as represented in the 
matrices X(t) and M{t), can be derived from the matrix exponential of an augmented 
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matrix. It can be shown 
matrix exponentials, 


that the matrix X(t) can be derived from a product of the following 


X(t) = e At [I 



-A B 

. 0 C. 




O' 

I. 


(48) 


Thus, computation of X(t) now involves the computation of a matrix exponential. A simple 
reliable algorithm for computing the matrix exponential is given in section 5. 

In a similar fashion, one can express the matrix M(t) in terms of a submatrix of a matrix 
exponential. To see this, we start from its definition 


M{t) = J j\ A[v ~ s) Be Cv De Es dsdv 

= / [J* e Av BeCv dt;} De Es ds 

-s:As;^} 


*De Es ds 


(49) 


Let’s perform a change of integration variable v — t r. We have, 

M(t) = S '} 

*De Es ds 

= - f e A(t ~ s) e~ Ar Be~ Cr dr^ 

*e Ct De- E{t ~ s) d(t. - s)e Et 

*e Ct De~ Eq dq e Et 


= / {[ e Mq - r) Be Ct/2 e- Cr dr j 

*e Ct/2 De E{t - q) dq (50) 


Notice that part of the integrand in equation (50) delimited by braces can be replaced by 
terms involving the exponential of an augmented matrix. This follows simply horn lesults 
developed lor the matrix X{t). With this substitution, we obtain 


M{t) = | [/ 0] exp 

%e Ct/2£) e E(t-<j) ^ 


A 

0 


Be Ct / 2 

-C 


0 

I 


(51) 
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= / V 0] exp 

JO 


A Be ct ' 2 
0 -C 


(t - q) 


0 

1 


e Ct l 2 D \ e Eq 


dq 


f 

' A Be Ct / 2 0 

} 


' 0 ' 

[I 0 Ojexp < 

0 -C e Ct > 2 D 


> 

0 

l 

0 0 E 

J 


I 


In this section, we have shown that the matrices X{t) and M(t) can be formulated in terms 
of the solutions of some matrix exponentials. Their evaluation depends therefore strongly 
on the accuracy and reliability of numerical methods for computing matrix exponential. We 
will present one such algorithm in section 5. However for computational expediency, special 
consideration must also be taken to ensure the efficiency of the overall scheme when the upper 
limit t is large and one of the matrices /l, C or D is unstable. Also one must economize 
memory requirements associated with high dimensionality of the augmented matrix when 
computing the matrix exponential. These considerations will be elaborated in sections 6 
and 7 where we give precise algorithms for the computation of the matrices X(t) and M{t) 
respectively. 


5, Numerical Method for the Matrix Exponential 


Several numerical methods are available for the computation of the matrix exponential [11]. 
Among all these, an approximation method based on Pade series is found to be satisfactory 
[12]. An important component in any numerical routine for matrix exponential is the scaling 
of the matrix argument prior to the series calculation. Due to the simple result that t Ai — 
(e^/ 2 ) 2 , a. scale factor in terms of powers of two (be 2 m ) is often used. In this scheme, one 
can recover the actual value of the original matrix exponential by performing m squarings 
on the matrix exponential of the scaled matrix. The index m is determined based on the 
desired size for the scaled matrix. In our algorithm, scaling is applied to the original matrix 
until its oo-norm ||A|| (X> falls below 1/2. 

As mentioned above, the preferred series approximation in our computation ol the matrix 
exponential is the Pade series. Let’s review some of the unique features associated with the 
Pade series for the case ol a scalar function On its most basic terms, it is a rational 

function of z of a preselected order that approximates the function J F(z). For a given choice 
of the order of the numerator (say N) and of the denominator (say M), the Taylor series 
representation of this Pade series must match the power series representation of T{z) for the 
first (N+M+l) terms. Namely, 


Hz) ~ p£(z) 




(52) 


In fact, the most common form of the Pade series is known as the diagonal sequence where 
the numerator and the denominator have the same order (i.e M = N). While it is known 
that the Pade series for the matrix exponential (i.e T{z) = e z ) converges only slightly faster 
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than the Taylor series for a scalar argument, the improvement is more significant tor matrix 
argument. In the matrix case, Pade series involves computation ol a numeratoi matiix 
A [{At) and of a denominator matrix V{At). For a diagonal Pade series of order N, we have 


V(At) 


+ 

+ 

+ 


t pN-l)lN± At 
(2Ny.(N-l)'. L 

(2A>-2pAn m ,v 2 . ... 

(2N)'.2\(N-2)'A' > ^ 

(2N-))! N \ 


(2N)\i<(N-iy. 

N ^(Atf 


(Aty + 


(2N) 1 . 


(53) 


and 


M(At) = 
+ 
+ 
+ 


T _ (2jV-l)!N! M 

1 (2lV)!(N-l)’ 

2N-2j l_ ^ {M) 2 _ 


(2N)!2!(N-2)! ' 

1—1 V ( 2 N-t)!N! 1 /jM 1 4- • • • 

(-i rjm( At ) N 


The matrix exponential is simply given by 


.At 


V~ l {At) N (At) 


(54) 


(55) 


Invertibility of V(At) is usually ensured by proper scaling of the matrix argument At. 

Another important consideration in the Pade series is its length A'. Assuming that the 
matrix At has been scaled such that H/U^ is less than 1/2, the parameter N can be choosen 
according to [12] such that 

< W! > 2 < £ (56) 


-.3-2 N 


(2N)\ (2N + 1)! 


where e is a given desired tolerance for accuracy. 

With a Pade series of N terms where N is determined from above, the approximation can 
be thought of as the exact calculation of a matrix exponential for a ” nearby” matrix ( At + E) 
where E is the error matrix with ||£||c* < c||Af||co. The relative error of the approximation 
is bounded by the following inequality, 


| |e (^E)- e *|| oo 


l^lloo 


< tWAtW^e 




(57) 


Thus, reducing the oo — norm of the matrix At would indeed improve the numeiical accuiacj 
of the matrix exponential. It has also been shown that methods by series appioximation yield 
better accuracy if the matrix argument has been preconditioned. Additional improvement 
may therefore be gained by first preconditioning the original matrix. Anothei immediate 
benefit of lowering the oo-norm of the matrix being exponentiated is that the actual scaling 
factor m needed would also be smaller; thereby resulting in a fewer number of matrix mul- 
tiplications in the squaring procedure. As usual, preconditioning a. matiix tends to biing 
singular values of that matrix closer together (i.e. lower the condition number), thus avoiding 
situation where scaling factor is predominantly determined by a few large singular values, 
and causing significant loss of precision related to the set of small singular values. The most 
common method used in the precondition of a matrix is the Osborne s method [14], which 
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minimizes the Frobenius norm of that matrix (and thus indirectly lowering its oo-norm). 
However, extensive tests conducted so far seems to indicate that preconditioning of a matrix 
did not yield significant reduction in the oo-norm and a smaller scaling factor to justify the 
added computational efforts incurred in the Osborne’s method. The procedure of precondi- 
tioning a matrix is nonetheless recommended from the point of view of improved accuracy 
(see [15] and[17] ). 

In the implementation of our design algorithm for optimal low-order controller synthesis 
[9], a value of e = 10 -8 has been selected requiring therefore a 4-terms Pade series (i.e N = 4) 
in the evaluation of the matrices X l (t), C l (t) and M l {t) of equations (39)-(41). Additional 
considerations in the implementation of the proposed method for computing X(t) and M ( t ) 
are given in sections 6 and 7 . 


6. Detailed Algorithm for Computing A(t) 

As seen in the previous section, the matrix X{t) can be evaluated in terms of a matrix 
exponential as shown in equation (48). Conceptually, it is a simple and straightforward 
procedure to compute the matrix exponential ol any arbitrary matrix using the Pade series 
discussed in section 5. However, it becomes a nontrivial task when we try to implement an 
efficient algorithm that examines carefully the issues related to accuracy, speed and memoiy 
requirement. The basic difficulties lie in the fact that the matrix exponential is for an 
augmented matrix of a particular form, 


/ 

-A B 

,\ 

e~ At e~ M X{t) } 

exp< 

0 C _ 

J 

o 

a 

o 


where A = C T = F' ! + o‘7 according to our problem in equations (39) and (40) for the 
matrices *'(*/) and C{t s ) respectively. Clearly, if the system matrix A is stable (i.e all 
the eigenvalues of the matrix A have negative real parts) then one could easily encountei 
numerical overflow when evaluating the term e~ At even though the matrix integrals A (<) and 
jC(t) of interest are perfectly well-behaved. The overflow problem occurs most likely in the 
final squaring process. To arrive at a feasible approach in the evaluation ol X{t), one needs 
to examine in details the steps taken in arriving at the matrix exponential of the original 
matrix starting from that ol a scaled matrix (i.e in the squaiing piocess). 

Let’s assume that one has scaled the input matrix A by AAt where At is a reasonably 
small time interval given by At = tjn = t/2 m . Thus, we need to first evaluate 


exp 


-A B 
0 C 



where At = t/n 


t/2 m . 


For notation convenience, we define 

D E ‘ 
0 F 

\ e ~ AAi t~ AAt f 0 At e Ar Be Cr dr 


exp 


-A B 
0 C 


At 


(59) 
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Furthermore, let W = exp{AAt) = D~ l . Now we can write our result as follows, 


X(t) = W n [D n ~'E + D n ~ 2 EF 

+D n ~ 3 EF 2 + h EF n ~ l ), 


or 


X(t) = W[E + WEF + W 2 EF 2 + • • • + W n l EF n ! ]. 

The above results are produced by performing m squarings of 


( 60 ) 


(61) 


D E 
0 F 


and taking the 
0)), the solution 


appropriate submatrix for X(t). In our application (cf. equations (39)-( 2 
would therefore involve products of matrices of size 2(n + r + n'). Close examination of 
equation (61) leads to the following algorithm involving only product of matrices of size 
_|_ r with the final result achievable in m steps, 


Step 1 : 

P\ = w. 

ii 

<5 

Rx = 

F 

Step 2 : 

Pi = Pi 

Q'2 — Qi + 
P\Q\R\, 

Ri = 

R 2 

Step m : 

p — pi 
1 m x m — 1 

Qm — Qm — 1 d - 
P m -lQm-1 Pm — 1 1 

Pm 

Rl-X 


Finally, X(t) = WQ m . It should be noted that one can "’absorb” this extra factor of W 
( = e AAt ) into the matrix Qi without any change to the above algorithm (i.e starting the 
above algorithm with Q, = WE instead). This removes the need to retain the matrix IT 
throughout the computation. 

Finally, one notes that the terms P t or R x for (i = l,m) may underflow and become a 
null matrix for some *; in particular when the scaling factor is large (i.e m large). When this 
situation happens, one can simply truncate the series calculation for X{t) up to the i th step 
in the above algorithm since all of the significant (and nonzero) terms have already been 
accumulated into the matrix Q{. 

7. Detailed Algorithm for Computing M ( t ) 

Here the numerical algorithm is a bit involved compared to the one given for the calculation 
of X(t). This is largely due to the increased complexity of the argument of the matrix 
exponential. Following the procedure described in section 6, let’s perform a scaling upon 
the input matrix A by A A t such that computation of the matrix quantities AA 0 , H • J, P, V 
and W = F _1 is well-behaved. These quantities are defined from the following matrix 
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exponential, 


f 

' A Be ct > 2 0 ] 


exp < 

0 -C e Ct ' 2 D 

At > 

O 

o 

tq 

j 


f P He ct ' 2 Mo 

1 


( 62 ) 

I 

0 V e Ct/2 J 

0 0 U 

Due to the possible numerical underflow in the matrix e c</2 for large t, the matrices H and 
J are computed directly from the following definitions, 

rAt 

/ 

JO 


c A 'U,: C ' dre CA ‘ 


(63) 


and 


J = e 


-CAt 


r 


e CT De KT dr 


E T 


(64) 


However, the computation of M 0 in equation (62) can still underflow due to its explicit 
dependence on e Ct H . For the calculation of the matrix M(t), ideally it can be obtained from 
m squarings of equation (62). If carried out in this manner, potential numerical overflow 
is eminent since, according to our equation for M l (tf) in (41), we have A —C — E — 
fi -|- a'L Hence, if the matrix C is stable, then the matrix exponential e = V n will 
become unbounded. To bypass this difficulty, as in the calculation for X{t), one needs to 
conduct the squaring algorithm explicitly. It can be shown that the matrix M(t) can :>e 
computed as 

P n -'M 0 + P n ~ 2 M 0 U 
P n ~ 3 M 0 U 2 + • ■ • + PM 0 U n ~ 2 
M 0 U n - 1 (65) 

HW 2 J + HW 3 JU + P HIV 3 J 
PHW 4 JU + P 2 HW 4 J + • • • 

P n ~ 2 HW n J 

This formulation no longer involves the matrix V. The above series lor M can be distin- 
guished into two parts— one that contains the matrix M a and the other that does not. The 
terms involving AA 0 can be thought of as 


M(t) 


+ 

+ 

+ 

+ 

+ 


[I 0 ] 


‘ P 

Mo ' 

n 

' 0 ' 

0 

U 


I 


( 66 ) 


which can be performed by m squarings. The remaining terms involving H, J. W, P, and U 
are computationally intensive and are of the form 


y P' HW 2+,+ i JU 2 , where 2 + i + j < 

i- 0 j=0 


n. 


(67 


This equation, owing to the restriction 2 + i+j < n, is not easily calculated in m(- log 2 {n)) 
steps. A reasonably efficient procedure for computing the final matrix M{t) is to merge both 
the easily computed portion given in equation (66) and the more difficult series in equation 
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(67) into a sequence of m steps, as shown in figure 3. Due to potential numerical underflow, 
the term W'~ 2 is not accurately obtained from the product W' V 2 where V = W~\ Indeed 
one needs to recompute the term W n+2 ~ v at each step of the above algorithm. This could 
become the major drawback of our scheme even though we have used an efficient matrix 
exponentiation routine to compute W x requiring at most 2 * log 2 {i) matrix multiplies. 

If in addition W n is zero (or effectively so), restriction on the indices i and j of 2+i+j < n 
in equation (67) becomes inconsequential; hence we can express 

E’i=o £ "=o pi H W 2+i+] JW = 

(. H + PHW + P 2 HW 2 + --)W 2 
*(J + WJU + W 2 JU 2 + • • •) 

resulting in a simpler algorithm involving the following three terms: 


■ p 

Mo ‘ 

n 

' 0 ‘ 

0 

U 


I 


( b ) (H + PHW + P 2 HW 2 + • • • + P n ~ 2 HW n ~ 2 ) 

(c) ( J + WJU + W 2 JU 2 + • • • + W n - 2 JU n ~ 2 ) 

This algorithm can again be computed in rn steps as seen in figure 4, but now there is no 

costly evaluation of a W k term at each step. 

Further simplification of the above algorithm can be achieved if we make use of the fact 
that we have A = E = C' T (cf. equation (41)) and therefore U = P = W T . If, lor some 
index j < m, W v (and likewise U v and P v ) is zero or nearly so, then this calculation for 
Ai{t) is reduced to .Vf ( t ) = HjW 2 Jj since M m - 0 , Hj = II,,, ■ and •/, — J m - 

In the following sections, we compare the usefulness of the proposed algorithm to the 
early algorithm presented in [9] in the design of low-order optimal controllers for two flexible 
mechanical systems . 

8. A Simple Two-Mass-Spring Design Problem 

Control of flexible mechanical systems has been of interest in recent years [18]. This problem 
provides us a simple design case where degeneracy in the closed-loop eigensystem can be 
easily illustrated. The problem is to control the displacement of the second mass by applying 
a force to the first mass as shown in Figure 5. At the start, it is simple to verify that the 
basic open-loop system has a pair of degenerate eigenvalues at the origin. Equations ioi the 
dynamic model are given below, 

miy\ = k(y 2 - yi) + u + w /ggj 

m 2 y 2 = - y 2 ) 
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or 



' Vi ' 


1 

O 

0 

O 

1— H 


' Vi ' 

± 

y[ 


0 -k/mi 0 k/mi 


y\ 

dt 

V2 


0010 


y 2 


. y'i . 


0 k/m 2 0 — k/m 2 


. y'7 . 


+ 


0 

1/mj 

0 

0 


( 69 ) 


(u -f iu) 


where m\ = m 2 = k\ = k 2 = 1* For comparison, we have obtained an optimal second-order 
controller design of the form, 


A = 


C 


0 1 

A21 A22 


B 


[Ci 1 c 


12 


D = [D 


11 


m 


using both algorithms. The control design problem is to minimize the following H 2 - nor 
of the closed-loop transfer function T yrW between the disturbance w and the displacement 
1/2 of the second mass through the controller design parameters A 21 , A 22 , Cn, U 12 and D\\. 
We start with the following arbitrary initial design guess of A 2 1 = -2,A 22 = ~UC n = 
0,Ci2 = 0.5, and D n = 0. Both algorithms converge effectively to the same optimal design 
gains of A 21 = -0.8571, A 22 = -0.9258, C n = 0,C 12 = -0.4535 and D n = -0.2449 and 
with an optimum value \\T y2W \\l = 7.71838215122. A summary of the resulting closed-loop 
eigenvalues is given in Table 1. The main difference between the two algorithms is in the 
CPU time for the overall computation. Results are obtained for a VAX/VMS- Workstation 
DEC-3500 as follows: CPU time of 19.59sec with the algorithm based on diagonalization 
and 97.36sec using the proposed method. The increase in computational load is expected 
and constitutes the basic trade-off between reliability and speed of the solution algorithm. 
The proposed algorithm is more reliable and with this advantage does take a. bit longei in 
computational time. With the early algorithm of [9], one cannot initiate the seaich foi an 
optimal compensator design with zero gains (i.e A 2 \ — A 22 — Cn = C\ 2 = D u — 0) because, 
in this case, the closed-loop system would have two pairs of degenerate eigenvalues at the 
origin; one for the rigid-body mode and the other from the open-loop compensator poles. 
To alleviate this problem, it was suggested that one simply starts with any compensator 
design (stabilizing or not) that produces initially a non- degenerate closed-loop system. Even 
with these considerations, it was found that occasionally the algorithm could break down 
due to the presence of near degeneracies in the closed-loop system matiix. 1 hus, loi a 
reliable design method, solution algorithm must treat degeneracies as a common occurrence. 
This situation is more evident in the optimal output-feedback control design lor high-order 
structural models with closely packed flexible modes. 

Future consideration would be to develop a hybrid algorithm taking advantage of the 
computational efficiency of diagonalization when the closed-loop system matrix is not de- 
generate, and turning to the current algorithm when degeneracies are detected. System 
degeneracy can be easily checked from the condition number of the eigenvector matrix. 
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In the next section, we present a design problem where degeneracies occur frequently and 
therefore it could pose a serious difficulty for the early design algorithm based on diagonal- 
ization. 


9. The JPL Large Space Structure Control Design 

In control problems for large flexible mechanical systems such as space structures, causes of 
eigenvalue degeneracies are usually more subtle in nature than the simple case presented in 
section 8 for a two-mass-spring system. The JPL large space structure has been caiefully 
designed to simulate a lightweight, non-rigid and lightly damped structure in a weightless 
environment [16]. The structure itself resembles a large antenna with a central boom-dish 
apparatus and an extended dish consisted of hoop wires and 12 ribs (Figure 6). There are 
two torque actuators (labelled H Al and HA10) on the boom and dish structure to control 
the two angular degrees of freedom in pointing maneuver, and force actuators at four rib 
root locations (labelled RA1, RA4, RA1 and flAlO) for vibration control. From the point 
of view of control design, it is a challenging problem since the plant has many closely spaced 
modes and is of reasonably high order. There are a total of 30 modes in the basic stiuctuial 
model. The flexible modes are lightly damped with damping ratios ranging from 0.007 to 
0.01. The two rigid-body modes have a damping ratio of 0.12. Our design concept is to 
use two available angular displacement sensors HS 1 and HS 10 of the boom-dish apparatus 
and the two torquers H Al and H A10 collocated with these sensors for control synthesis. 
With this selection, 20 of the flexible modes associated primarily with the rib motion become 
uncontrollable and unobservable. These modes are removed by modal truncation fiom out 
plant synthesis model. Eigenvalues of the remaining 10 modes are shown in Table 2. 

An optimal low-order controller is designed to dampen vibration ol the antenna to exter- 
nal excitations. To evaluate the effectivenes of the control system, we perform the following 
test. The entire structure is agitated using the two boom-dish actuators for the first 6.4 sec- 
onds with an applied torque in the form of a square wave of 0.8 second in width and with an 
amplitude of 1 N-m. The control system is then activated right after the excitation has been 
removed, and responses of the excited structure at the sensors are examined. The design 
objective is to damp out the induced vibration as fast as possible without excessive use of 
controls. Note that the natural responses of the structure will take about a few minutes to 
decay to zero (Figure 8). 

For practical implementation, the controller design is choosen to be of 6* order and has 
the following form, 

—50 0 Ai3 A\4 A 15 

0 —50 A 2 3 A 24 A 2 5 

0 0 0 1 0 

0 0 A43 A44 0 

0 0 A53 A 54 0 

0 0 Te3 .4(34 /If; 5 


A 16 
A 26 
0 
0 
1 

A 66 


70) 
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B 


C 


B u 

B 12 

B 21 

B 22 

B 3 I 

B 32 

B 4 I 

B 42 

B$\ 

B 52 

Bq\ 

Bg 2 . 


50 0 0 0 0 0 
0 50 0 0 0 0 


The first two states in the controller model serve as roll-filters, limiting the control bandwidth 
to less than 50 rad/ sec. In the design optimization, we have a total of 28 design variables: 
16 in the controller A matrix and 12 in the B matrix. The objective function for design 
optimization consists of a sum of weighted H 2 - norms of physical response variables observed 
at different location ol the structure. It is of the form 


J[tf) = 12 

Lim I ( £ QiE a [!/?((/)] + E R,E„ [«](*/)] 

t f -^CO Z [p = l j = 1 J 


Note that the expectation operator E a [~] is for a system destabilized by a factor o. Table 
3 lists the design variables y t and their corresponding penalty weightings Q % . Also given in 
the table are the control design weightings R, for the actuators H A\ and HA] 0. Responses 
in the above objective function are evaluated to random disturbances of unit white-noise 


spectra applied simultaneously at all the hub and lib actuatois. 

The design optimization begins with the following arbitrary initial guess on the controller 

matrices A and B , 


A 


-50 

0 

0 

0 

0 

0 


0 

-50 

0 

0 

0 

0 


1 

0 

0 

-2 

0 

0 


0 

0 

1 

-1 

0 

0 


0 

1 

0 

0 

0 

-4 


0 

0 

0 

0 

1 

-4 


B = 


0.1 

0 

0 

0 

0 

1 


0 

0.1 

0 

1 

0 

0 


A destabilization factor a ol 0.071 was used to ensure that all the closed-loop eigenvalues 
have a real part less than -0.071. The optimization fails to converge when a destabilization 
factor of greater than 0.075 was selected. This difficulty seems to be in moving the modes 
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at 1.68 Hz under this controller configuration, implying that additional degrees of freedom 
must be added to the controller structure given in equation (70). 

While the optimization convergence itself took 13.5 hours on a VAX/VMS Workstation 
DEC-3500, the proposed algorithm for the calculation of the objective function and its 
gradients with respect to the design parameters is robust and leads to well-behaved design 
convergence. The final optimal values of the A and B matrices are shown in Figure 7. 
Closed-loop eigenvalues are given in Table 4. Primary improvement is seen in the increased 
damping of two modes at 0.65 Hz. 

Closed-loop responses of the sensor and control variables corresponding to this design are 
shown in Figure 8. The controlled responses decay to zero in about 20sec after the excitation 
has been removed. Notice that the control torques are within the desired limits of 1 N-m; 
the results are obtained through adjustment of the control design weights Rj in Table 3. 
This design example demonstrates the usefulness of a design algorithm for robust low-oidei 
controllers using parameter optimization, and the accompanying improvement of solution 
reliability using the algorithms described in sections 6 and 7 for degenerate systems. 


10. Conclusions 

Numerical algorithms for computing matrix exponentials and integrals of matrix exponen- 
tials have been developed to handle cases where the system matrix is degenerate. Numerical 
optimization combined with the given algorithms for the evaluation ol the cost function 
and its gradients with respect to the controller design parameters has well-behaved conver- 
gence even when the closed-loop system becomes degenerate. These algorithms have been 
incorporated into a computer-aided-design package for synthesizing optimal output-feedback 
controllers. Reliability of the algorithm has been demonstrated using typical design prob- 
lems encountered in the control of flexible structures. Clearly this algorithm when combined 
with a previous one based on diagonalization would enhance significantly the overall reliabil- 
ity of the optimal design procedure for low-order controllers, thereby providing an effective 
automated design environment for multivariable control synthesis. 
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Figure 3: An m-Step Calculation of M(t) 
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Figure 4: A Simplified m-Step Calculation of M{t) 





Figure 5: A Two-Mass-Spring Mass System 


Eigenvalue 

Damping 

Freq (Hz) 

-0.2290 ± 0.3397* 

0.559 

0.0652 

-0.1553 ± 0.8480* 

0.180 

0.1372 

-0.0786 ± 1.2950* 

0.061 

0.2065 


Table 1: Closed-Loop Eigenvalues of the Two-Mass-Spring System 


Eigenvalue 

Damping 

Freq (Hz) 

-0.09500 

± 

0.7860* 

0.120 

0.12600 

-0.08575 

± 

0.7093* 

0.120 

0.13704 

-0.02802 

± 

4.0024*' 

0.007 

0.63701 

-0.02929 

± 

4.1844* 

0.007 

0.66598 

-0.07405 

± 

10.583* 

0.007 

1.68434 

-0.07405 

± 

10.583* 

0.007 

1.68434 

-0.11310 

± 

10.616*' 

0.007 

2.57123 

-0.11785 

± 

16.384* 

0.007 

2.67929 

-0.21365 

± 

30.520*' 

0.007 

4.85749 

-0.21365 

± 

30.520* 

0.007 

4.85749 


Table 2: Open-Loop Modes of the Antenna. Structure 
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Variable 

Q, 

Description 

mu 


Rib #1 root velocity 


3950 

Rib ^4 root velocity 

RSI 


Rib #7 root velocity 

RS 10 


Rib #10 root velocity 



Hub angular velocity 


ItMil 

Hub angular velocity 

■ 

1100 

Rib #1 root displacement 

RS4 

1050 

Rib #4 root displacement 

RS7 

1150 

Rib #7 root displacement 

RS10 

1025 

Rib # 10 root displacement 

HSl 

3900 

Hub angular displacement 

HS10 

4100 

Hub angular displacement 

Variable 

Ri 

Description 

HAl 

41 

Hub torque actuator 

HA10 

40 

Hub torque actuator 


Table 3: Design Variables for JPL Antenna Structure 


Eigenvalue 

Damping 

Freq (Hz) 

-0.086899 

± 

0.6588i 

0.1308 

0.1058 

-0.089071 

± 

0.741 Oi 

0.1193 

0.1188 

-0.3165 

± 

3.624i 

0.0870 

0.5790 

-0.2528 

± 

3.790i 

0.0666 

0.6045 

-0.2162 

± 

4.1 12i 

0.0525 

0.6553 

-0.2056 

± 

4.185S 

0.0491 

0.6669 

-0.074193 

± 

10.58i 

0.0070 

1.684 

-0.074589 

± 

10.58i 

0.0070 

1.684 

-0.1168 

± 

16.15i 

0.0072 

2.570 

-0.1253 

± 

16.83i 

0.0074 

2.678 

-0.2142 

± 

30.52i 

0.0070 

4.857 

-0.2143 

dt 

30.52i 

0.0070 

4.857 

-49.99 



1.000 

7.956 

-49.99 



1.000 

7.956 


Table 4: Boom- Dish- Controller Closed-Loop Modes 
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Figure 6: Antenna Structure 
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HA1 Control (N-m) HS1 Response (rad) 


0 2.874 

-2.270 

1.7322 x 10- 3 

50 1.225 

0.7825 

6.551 

0 0 

1 

0 

0 -15.73 

-0.8799 

0 

0 1.560 

0.2256 

0 

0 2.400 

-1.269 

-13.62 

r 

5.343 

-1.2310 x 10~ 4 

, 6.2118 x 10“ 4 

4.783 


2.701 

-8.1595 x lO" 4 

B = 

2.221 

9.3152 x 10" 4 


-0.5147 

5.379 


1.614 

-1.252 


-2.2131 x KT 4 
-1.037 
0 
0 
1 

-0.9810 


Figure 7: Optimized Controller Matrices for LSCL Problem 
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Figure 8: Open-Loop (solid curve) versus Closed-Loop (dashed) Responses 








DESIGN CHALLENGES FOR THE UH-60 ROTOCRAFT 


Brett VanSteenwyk and Uy-Loi Ly * 
Department of Aeronautics and Astronautics, FS-10 
University of Washington 
Seattle, WA 98195 


Abstract 

High performance rotocraft controller design is characterized by having to compensate both 
lateral and longitudinal dynamics in a single controller design. Rather than have these modes 
largely decoupled in the natural state, often one needs to incorporate mode decoupling as a 
part of the controller design, especially with the modelling of higher frequency dynamics. There 
exists the usual concerns of stabilizing the behavior without excessive actuator output, however, 
in a high performance rotocraft such as the UH-60, good response bandwidth is as important 
as stabilization and decoupling. 


1 Introduction 

The challenge is to develop a stabilizing output feedback controller using a 23 state model to 
simulate the UH-60 rotocraft in the hover flight condition. The essence of this challenge is to 
develop with a consistent set of techniques that allow a design to meet or surpass flight criteria. In 
addition to a stabilizing and decoupling controller, command following on <j> (pitch), 9 (roll), r (yaw 
rate), and w (vertical rate) is required. Sensor outputs include these terms, along with pitch and 
roll rates, and the forward and transverse velocities (8 sensor inputs total). Optionally, one could 
control the yaw angle (heading) instead of the rate — this issue is not a fundamental one insofar as a 
design approach is concerned. Techniques for designing controllers operating with heading should 
not be distinguished from those operating on yaw rate. 

The UH-60 rotocraft open loop model in the hover flight condition is mildly unstable (see 
Table 1). The pair of unstable poles represent a phugoid-like response in the front/side velocity 
coupling in with the pitch/roll. In addition, there are modes that are either pure integrating (yaw 
from yaw rate, for instance), or are near integrating (roll). The sensor output of the model consists 
of the angular rates p, q, and r, corresponding to the angles 9 (pitch), <p (roll), and V ; (yaw) which 
are also considered part of the sensor output. In addition, the translational velocities u, v, and w 
are considered sensed. 

There are 4 actuator inputs: <5o, <5 S , <5 C , and These stand for the collective (main rotor), 

the main rotor sine, main rotor cosine, and tail rotor pitch inputs, respectively. The response of the 
UH-60 to a unit collective pulse is shown in figure 1. These graphs give one a sense of the scaling 
and sensitivity of the various inputs and output terms. 

*The work of B. VanSteenwyk and U. Ly is supported in part by NASA Ames Research Center under grant 
contract NAG-2-691. 
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Eigenvalues 

Damping 

Freq (rad/sec) 

-9.38957 

± 

51.2473i 

0.180 

52.1003 

-5.74275 

± 

37 . 1461 i 

0.153 

37.5873 

-8.62407 

± 

24.48631 

0.332 

25.9606 

-24.5088 

± 

2.903361 

0.993 

24.6801 

-4.22353 

± 

19.6367i 

0.210 

20.0858 

-19.3508 



1.000 

19.3508 

-5.71403 

± 

5.43266i 

0.725 

7.88441 

-3.26916 

± 

4. 262801 

0.609 

5.37205 

-4.82826 



1.000 

4.82826 

-1.35460 



1.000 

1.35460 

0.229236 

± 

0.428933i 

-0.471 

0.486346 

-0.135879 

± 

0.526819i 

0.250 

0.54406 

-0.230418 

± 

0.01 lOOOi 

0.999 

0.23068 

0.00000 



0.000 

0.00000 


Table 1: Open Loop Eigenvalues of UH-60 Rotocraft 


2 Basic LQ Design 

What are the limits of full state feedback? Although full state LQ is not always the most satisfactory 
design technique, in an ideal circumstance it can define the limits of a stabilizing output feedback 
controller. Suppose one defines the control penalty matrix as the identity and focuses on the most 
effective output variables to penalize. In addition, there is the issue of controlling yaw rate, or 
the yaw (heading) itself. Unless noted otherwise, yaw rate is the controlled quantity (this will not 
change the basic conclusions). 

We initially try blindly to penalize all outputs to try to stabilize the system. Given a stabilized 
system, a feedforward matrix is constructed to try to input commands. It quickly becomes clear 
that the quality of command response to pitch and roll commands drops if one is trying to stabilize 
the u and v responses. This seems physically reasonable (pitch forward to increase forward velocity, 
and so on). Given an R of I, it seems that increasing the response weighting Q will not change 
the least stable eigenvalues much-they go from —0.71278 ± 0.2257? to —0.71289 ± 0.2251?' when Q 
is raised from 10*1 to 1000*1. A full listing of the eigenvalues of the closed loop system is given 
in Table 2 for the case where Q = 10*I. Although the slowest poles are not that slow, one has an 
essentially degenerate pair, and thus one has their persistent behavior. It is unfortunate that they 
do not move with increasing Q. 

It turns out that one can control pitch and roll directly, and thus the ground velocities u and 
v as a consequence, or one can control the ground velocities with pitch and roll control derived 
from this. Since pitch and roll are more desired as directly commanded, and it is more reasonable 
to use these to control ground velocities. The penalty weighting on u and v is dropped. The 
resulting optimization produces a set of stabilized plant eigenvalues that have worse modes than 
before; however, it seems that these modes stay much less disturbable and observable from the 
commanded states. Increasing Q reduces the interaction with these poles without changing their 
values-they become more exclusively associated with the u and v states, and less with the other 
responses. Even so, as seen in the command responses, there is some '’hangofP even when the 
nominal diagonal entries for Q (except u and v) are 100,000 (see Table 3). Another way of viewing 
this hangoff is to observe Bode plots of pitch relative to pitch command, etc. (See figures 3 and 4) 
With the zero frequency dropoff in amplitude of response, one would see the dynamics manifested 
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Figure 2: Basic LQ Full State Design Approach 


Eigenvalues 

Damping 

Freq (rad/sec) 

-33.3483 

± 

18.5324 i 

0.566 

58.8855 

-10.9609 

± 

53.79Ui 

0.200 

54.8964 

-41.2844 



1.000 

41.2844 

-12.3645 

± 

28.4 1 1 7i 

0.399 

30.9856 

-2.97880 

± 

17 .3773i 

0.169 

17.6308 

-29.3151 



1.000 

29.3151 

-23.0359 

± 

4.57704 

0.981 

23.4862 

-18.4648 



1.000 

18.4648 

-9.25176 

± 

5.21430 

0.871 

10.6200 

-11.4372 



1.000 

11.4372 

-9.81321 



1.000 

9.81321 

-5.15694 

± 

3.57155i 

0.822 

6.27295 

-0.72756 

± 

0.17 456i 

0.972 

0.74821 

-0.71278 

± 

0.22572i 

0.953 

0.74767 


Table 2: Closed Loop Eigenvalues of UH-60 Rotocraft, All Output Penalty, Q — 10*1 





Phase (Degrees) Magnitude (dB) 


Poles 

Damping 

Freq (rad/sec) 

-4558.04 

1.000 

4558.04 

-1537.43 

1.000 

1537.43 

-560.935 

1.000 

560.935 

-104.625 ± 88.1763i 

0.765 

136.826 

-89.8398 

1.000 

89.8398 

-74.8621 ± 22.1033i 

0.959 

78.0570 

-11.5621 ± 53 -4352i 

0.211 

54.6718 

-29.4031 

1.000 

29.4031 

-28.8883 

1.000 

28.8883 

-2.57737 ± 17. 1876i 

0.148 

17.3798 

-10.6260 

1.000 

10.6260 

-10.2186 ± 0.749891i 

0.997 

10.2461 

-4.75180 ± 5.70887i 

0.640 

7.42770 

-1.03407 

1.000 

1.03407 

-1.00059 

1.000 

1.00059 

-0.965363 

1.000 

0.965363 

-3.98898e-003 ± 6.71659e-003i 

0.511 

7.81182e-003 


Table 3: Poles in LQ Stabilized UH-60 Rotocraft System, Q= 100,000 



Frequency (Rad/Sec) 



Frequency (Rad/Sec) 


Figure 3: Bode Plot of Pitch/Pitch Command Response, LQ for Q— 100,000 
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Figure 4: Bode Plot of Roll/Roll Command Response, LQ for Q_ 100 ,000 

in a hangoff. Note that this problem is not very severe for the yaw and up rate commands, though 
noticeable. 

A ready conclusion was formed-something other than straight penalty on the command outputs 
was needed. To get decent results, the value of Q (and subsequent actuator activity) had to be 
stratospheric. There were still minor, but noticeable, hangoffs in the commanded outputs even at 
values of 100000. 

Note the following meaning with u and v being nearly integral poles. It simply provides for 
their control via pitch and roll commands. To attain a desired forward velocity u, one would pitch 
forward, allowing the rate to ramp up to a desired value, then pitch level again. It would seem 
intuitive, then, that it is impossible to regulate the set: [u,v,<p,0] as a whole. It would also seem 
that the existance of near- integrator poles for u and v is desirable. 


3 LQ Design with Integral Control on Pitch and Roll 

A ready solution to relieve the hangoff in a command response is to put an integral control state 
on it. Although for the purpose of synthesis we put the integrators in the plant, m reality they are 
more states in the controller. Thus any practical implementation will have to contend with these 
additional states in the controller and the additional numerical challenges they provide. Addition- 
ally, integral control is associated with lower bandwith, and thus is to be avoided if possible. It 
seemed that the pitch and roll commands had the worse of the hangoff problems, so integral control 
was applied for these commands. In addition, in the interest of eliminating another integrator, the 
heading state was not included (that the pilot controls heading via yaw rate commands). 

The resulting system had 25 states: 23 model states, and 2 integrator states. The existance of 
integral states allowed the blending of proportional, integral, and derivative states for both pitch 
and roll to form zero locations that would serve as closed loop system pole attractors. In a frequency 






Figure 5: LQ Full State Design with Two Integral Controls 
domain sense, one had: 

a\ * J N) + ci2 * A0 4- as * 9 — - * ^ cti + c 12 $ + n 3 s 2 ) 

To keep the scaling relative to the error constant, one should leave a 2 (the proportional term) as 1. 
Thus, ai and a 3 need to be adjusted so that the zeros associated with this criterion (penalty) will 
be at desired locations. For the pitch criterion, — 1 , a 2 = 1, and a 3 = 0.25. For the roll criterion, 
ai = 1, a 2 = 1, and a 3 = 0.2. For Q - diag[l ,0.1 , 1.6,2] ( 9 , (f>, r, and w), and R - diag[10, 1,1,1] 
(<5o, 6 S , b c , and Str), one achieved very satisfactory results, revealed in part by the behavior of the 
eigenvalues (see table 4). The Bode and time series responses tell the rest of the story (figures 6 to 
13). 

4 Loop Transfer Recovery 

Given the selection of LQ design as the basic approach for the full state case, the use of LTR 
to select appropriate observer dynamics (full or reduced order) largely presumes an existing set 
of feedback gains. The degrees of freedom reside with the estimator portion (be it full state or 
otherwise). The preferred norms and criteria for recovering the full state performance from an 
output feedback controller should be kept consistent. Thus, since we have started by using we 
shall continue to do so. 
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Phase (Degrees) Magnitude (dB) 


Poles 

Damping 

Freq (rad/sec) 

-9.38043 

± 

51.2638i 

0.180 

52.1149 

-6.11002 

± 

37.2574i 

0.162 

37.7550 

-8.72514 

± 

24.601 li 

0.334 

26.1025 

-24.5081 

± 

2.91 104i 

0.993 

24.6804 

-4.20761 

± 

19 53691 

0.211 

19.9849 

-19.3829 



1.000 

19.3829 

-5.98385 

± 

5.24203i 

0.752 

7.95521 

-3.75254 

± 

4.674 19i 

0.626 

5.99413 

-6.51433 



1.000 

6.51433 

-4.65048 



1.000 

4.65048 

-1.87052 

± 

1 .379951 

0.805 

2.32446 

-1.48464 

± 

0.6005 1 i 

0.927 

1.60149 

-1.70573 

± 

0 . 1 0905i 

0.998 

1.70922 

-3.98898e-003 

± 

6.71659e-003i 

0.511 

7.81182e-003 


Table 4: Poles for 2 Integral LQ Stabilized System 
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Figure 6: Bode Plot of Pitch/Pitch Command Response, LQ 2 Integral Control 
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Phase (Degrees) Magnitude (dB) Phase (Degrees) Magnitude (dB) 
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Figure 7: Bode Plot of Roll/Roll Command Response. LQ 2 Integral Control 
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Figure 8: Bode Plot of Yaw /Yaw Rate Command Response, LQ 2 Integral Control 
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Frequency (Rad/Sec) 

Figure 9: Bode Plot of Up/Up Rate Command Response, LQ 2 Integral Control 
The premise of CLTR is simple. Assume the system: 


Ax 

+ 

B\w 

+ 

B 2 u 

(dynamics 

C\x 

+ 

D n w 

+ 

D 12 u 

(criteria) 

C 2 x 

+ 

D 2 \w 

+ 

D 22 U 

(sensors) 


where w is the command/disturbance input and u is the control input. A state feedback matrix 
F is found that gives the system desired properties. Since state feedback represents only the 
achievable, but not the actual performance due to noise and limited sensor output, the focus is to 
add dynamics to synthesize individual or combinations of the necessary states. Within the context 
of this problem, the disturbance input w is the commanded output (pitch roll, yaw, or vertical 
velocity), with u remaining as the full set of actuator controls. The objective is to minimize the 
difference in criteria responses to disturbance (command) inputs for the full state veisus output 
feedback systems (see [4]): 

T zw (s) - Tq zw (s) = Ej(s) = T zu {s)Mj(s) 

This is generally approached by minimizing ||M/(s)|| 2 in Ricatti-based methods, while the error 
function itself can be minimized using direct optimization schemes. 

In addition, the direct optimization does allow one to alter some of the state feedback gains. 
In this case, it was initially held that those states associated with <p and 0 and their error integrals 
are held constant as these nominally would not be affected by changes in the rest of the controller. 
However, for the states associated with the rest of the outputs (p, g, r, u, u, and w), the gams 
of the feedback matrix are allowed to float. Although the integrators on the command error are 
treated as part of the plant here, it is understood that they will be a part of the controller. The 
formal numerical approach induces one to leave it as part of the plant during the design process so 
as to use the full state feedback matrix as is. 
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Roll rate u r: Yaw rate Tail Rotor Roll Collective 
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Figure 10: LQ 2 Integral Control 
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: Roll rate u r: Yaw rate Tail Rotor Roll Collective 
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Figure 11: LQ 2 Integral Control 
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Roll rate u r: Yaw rate Tail Rotor Roll Collective 
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Figure 12: LQ 2 Integral Control 
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Figure 13: LQ 2 Integral Control 
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5 Output Feedback Designs from CLTR: The 2 Integral State 
System 

These arise by starting with the full state feedback matrix (or relevant parts thereof) and moving 
to reduce the error difference between the trial output feedback system and the baseline design 
produced in Section 3. What results can one expect? Looking at the S C B. [3] of the basic 
plant from the control input to the measured output, we find that it is left invertible and include 
4 infinite zeros. There are no invariant zeros. The infinite zeros are probably from the integral 
control states and combinations of u, v, pitch, and roll. Note that the original state feedback matrix 
Iv is preferable for constructing a command feedforward for r and w. While exact recovery is not 
possible due to some infinite zeros, asymptotic recovery is. The results from the full order and 
reduced order observer designs are expected to be good. Not much else offhand can be said for the 
general compensators. 

5.1 LTR Design I: Full Order Observer Design 

In regards to reducing this error, 

TUs)Mf(s) = Ej(s) = 7 l w (s) — Tq zu ,(s ), 

one can compute the error Mj(s ) directly in terms of the plant matrices and the estimator matrices. 
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Plant Structure 


Ob server Structure 

= (.4 - KC 2 )x + B 2 u + A'?/ 

= -Fx 


x = Ax + B\w + B 211 
z = C\x + D\\w + Duu 

y = C2X + £> 21 ™ + D22II 

With K is the unknown estimator gain matrix 

Mj(s) = F(«7 - .4 + I<C 2 r\B 1 - KD 2X ) 


M 


/ 


s) = A T + Cf/f T ) 

This corresponds to a dual system: 

X = T 7 x + Cju + F T ™ 

2 = where: 

u = -K t x 

Minimizing 2 above in an H 2 sense with respect to K T is most easily performed through LQ. This 
is in general a problem with coupled weighting between the states * and controls u, and is usually 

singular. 


7773 

A 


in f 00 r min r x r n 
f / z 1 z clt - A / [x ti] 

Jo Jo 



B\ * Z?2i 




£>21 * 

D 21 * D\i 


i 

1 


dt 


The issue is to deal with the singular part of D^\ ■ If one does the Singulai \ alue Decomposition 

OtI D 21 • /p r tT 

D 2 i — t/. d^d Vd ; A>i * £>21 = Ud^d^d^d 

where: 


'Zd'Z'd = 


r -i 2 

0 


0 

0 

0 


0 

<^2 


0 

0 


One must supplement the singular part of this matrix with a small magnitude matrix. The singular 
part, in practice, is usually not well defined: the singularity is usually relative to the maximum 
singular value, not just an element being zero. Suppose one chooses an c equal to 0.000001 *(T 1 . If 
this magnitude is greater than a 2 m for some m, then the matrix to add to £>21 * £>21 would be: 


0.000001 * cr\ * Ud 


[" [0] (m-l)x (m — 1) 

0 - 

1 

O 

Mi 


V? 


One examines the asymptotic properties as e — 0. If the open loop system can recover the state 
feedback properties exactly, then the gain matrix K 1 will converge. If not, usually some of these 
gains will approach infinity. 

For this design, K will not converge as e goes to 0. What was done here was to observe the 
recovery error as a function of e and figure a tradeoff between decrease in error and increase in 
controller activity. It became simpler than that when it was observed that the response to the yaw 
rate and vertical rate commands became worse at a point. This seemed to correspond to where 
the recovery error no longer improved a great deal with decreasing e (the responses to the pitch 
and roll commands were still improving at this point). A value at 0.0025 was settled upon. The 
following shows the recovery error lor some values of e. 
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e 

Recovery Error 

0.01 ^ 

0.58031 

0.005 

0.48800 

0.0025 

0.44232 

0.00125 

0.43586 


The resulting eigenvalues for the chosen design are in Table 5. These tend not to be very sensitive 
to the value of e chosen and hence do not tell very much. What does tell a lot is a plot of the 


Poles 

Damping 

Freq (rad/sec) 

-251.746 n 

1.000 

251.746 

-33.5097 ± 49 -0920i 

0.564 

59.4384 

-10.1434 ± 52.61431 

0.189 

53.5832 

-9.38043 i 51/2638i 

0.180 

52.1149 

-46.5182 

1.000 

46.5182 

-6.11002 ± 37.25741 

0.162 

37.7550 

-10.4599 ± 27.13901 

0.360 

29.0849 

-8.72514 ± 24.60111 

0.334 

26.1025 

-5.92371 ± 19.6943i 

0.288 

20.5659 

-4.20761 ± 19.5369i 

0.211 

19.9849 

-24.4609 ± 3.59786i 

0.989 

24.7241 

-24.5081 ± 2.91 104i 

0.993 

24.6804 

-19.3829 

1.000 

19.3829 

-20.0000 

1.000 

20.0000 

-20.0000 

1.000 

20.0000 

-11.1077 

1.000 

11.1077 

-5.98385 ± 5.242031 

0.752 

7.95521 

-3.84041 ± 5.56737 i 

0.568 

6.76345 

-3 75254 ± 4.67 4 1 9i 

0.626 

5.99413 

-6.51433 

1.000 

6.51433 

-3.56047 ± 3.39555i 

0.724 

4.92003 

-4.70424 

1.000 

4.70424 

-4.65048 

1.000 

4.65048 

-1.87052 ± 1 .379951 

0.805 

2.32446 

-1.48464 i 0.6005 li 

0.927 

1.60149 

-1.70573 ± 0.1 0905i 

0.998 

1.70922 

-016719 ± 0.49056i 

0.323 

0.51827 

-0.23179 ± 0.41 152i 

0.491 

0.47231 

-1.28631 

1.000 

1.28631 

-4.26079e-3 ± 6.65864e-3i 

0.539 

7.90518e-3 


Table 5: Full Order Observer Design Eigenvalues, 2 Integral Control 

maximum singular value of both the controller activity and of the error with respect to frequency. 
These are shown for the various values of e in figures 15 and 16. 

The design does not end here. Note that this recovery process operates on a matrix M, 
not the actual recovery error. Perhaps much of the controller activity has gone into regulating 
high frequency components of the error that will not pass through the plant anyway (remember: 
Ef(s) = Tj u (s)M f(s)) Suppose one can augment the transposed system associated with M to 
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Results for Various Weighting Factors 



Results for Various Weighting Factors 
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Figure 16: Controller Activity, Full Order Observer 





reduce weighting of the subsequent LQ process on the higher frequency components. Thus, the 
recovery design will do better at the low frequency end at the expense of the high frequency end. 
Intuition would indicate that putting the corners of these first order filters on the outputs at 8 ra- 
dians/s for the pitch and roll components and 4 radians/s on the yaw rate and up rate components 
would make sense. Actually, some testing was done, and a corner of 0.5 radians/s was used for all 
filters. Because of this, weighting for more controller activity (lower values of e) were calle or. 
Some overall results can be summarized as follows: 


6 

Recovery Error 

0.0001 

0.3629 

0.00001 

0.2777 

0.000001 

0.2746 

0.0000001 

0.4136 


Poles 

Damping 

Freq (rad/sec) 

-208.040 



1.000 

208.040 

-114.649 



1.000 

114.649 

-47.6094 

± 

52.9427i 

0.669 

71.2010 

-10.3964 

± 

52.4826i 

0.194 

53.5025 

-9.38043 

± 

51.2638i 

0.180 

52.1149 

-62.2535 



1.000 

62.2535 

-62.2535 



1.000 

62.2535 

-6.11002 

± 

37.257 4i 

0.162 

37.7550 

-13.5964 

± 

27 .2780i 

0.446 

30.4788 

-8.72514 

± 

24.601 li 

0.334 

26.1025 

-4.20761 

± 

19.53691 

0.211 

19.9849 

-7.14504 

± 

18.18571 

0.366 

19.5390 

-24.5250 

± 

3 77065i 

0.988 

24.8131 

-24.5081 

± 

2.91 104i 

0.993 

24.6804 

-19.3829 



1.000 

19.3829 

-12.0730 



1.000 

12.0730 

-5.98385 

± 

5.24203i 

0.752 

7.95521 

-3 19825 

± 

6.34656i 

0.450 

7.10687 

-3.75254 

± 

4.67419i 

0.626 

5.99413 

-6.51433 



1.000 

6.51433 

-3.96063 

± 

2.27644i 

0.867 

4.56823 

-4.88859 



1.000 

4.88859 

-4.65048 



1.000 

4.65048 

-1.87052 

± 

1.37995i 

0.805 

2.32446 

-1.48464 

± 

0.60051i 

0.927 

1.60149 

-1.70573 

± 

0.10905i 

0.998 

1.70922 

-0.66682 

± 

0.399321 

0.858 

0.77724 

-0.23421 

± 

0.464261 

0.450 

0.51999 

-1.39861 



1.000 

1.39861 

-4.26079e-3 

± 

6.65864e-3i 

0.539 

7.90518e-3 


Table 6: Full Order Observer Recovery, 2 Integral Control, Weighted 
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As one can see, this design is much more favorable. The errors at higher frequencies for M } are not 
significant in the overall error. Note that now there is a definite minimum at 0.000001 (the chosen 
e). The eigenvalues are seen in table 6, with plots of the recovery error and controller activity as 

seen in figures 17 and 18. 



5.2 LTR Design II: Reduced Order Observer Design 

If one were to reconstruct all the states assuming that the measurements were relatively noise-free, 
then a Luenberger Observer with a minimal number of estimated states would be appropriate. The 
variables: p, q, r, u, v, w, <j>, and 6 are all assumed to be measured (along with the integrals of <f> 
and 9), and the remaining 15 states, those of inflow, lag and flap, had to be estimated or otherwise 
accounted for. Thus, in the general form for the Luenberger Observer: 


v — Lv -f G\ y 4- G 2 11 
x ~ Pv 4 Jy 
u = -Fx 


with: 

QA — LQ = G 1 C 21 O 2 — QB 2 , and: JC 2 4- PQ — ^ 23 - 

A clever set of transformations exists that will move this problem into the form of a full order 
observer and a corresponding dual system (see [4], section 4.2). Fortunately for this problem, 
little has to be done to achieve this form. The direct feedthrough term from the controls to the 
sensor output is assumed zero. The first step involves a transformation of the original system. 
The measured outputs become states (unless there is a direct feedthrough from the disturbances to 
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Results for Various Weighting Factors 
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Figure 18: Controller Activity, Shaped Full Order 


them). The rest of the states along with the direct feedthrough terms form the rest of the system 
Original System: 

• i 

x = 

y - 

becomes: 

_ 
i ’2 

yo _ 

. XJi 

The UH60 system has neither a D 2 1 nor a Z? 22 term, so the x 2 states are soley represented by the 
non-output states. Needless to say, C’ 2 , 02 is non-existant (as is y 0 ), and no transformation on the 
state variables is needed since the existing outputs are states already. Defining an observer gain 
matrix K r , with partitions [A' r o, K r l] corresponding to [y 0 , S/i], » few more transformations will 
yield the Luenberger form: 

First define: 

A or — A22 — ^0^2,02 ~ kr\A \2 
Then: 

V — A or v + ( £$2,2 “ A rl ^21 ) U + A ,0 £/0 

+ (A21 — A r i An + A or A r i) yi 


A'x' + B[w + B! 2 u 

C\x + D}iW + Dj2^ 

CW + D' 


21 ‘ 


An 

A 21 


A\2 

A 22 


Xi 

X2 


0 


1 


p—m 0 


C 1 2,02 

0 


+ 


Ci a; + D\\w + D\ 2 U 


x x 

X2 


B 1A 

B\ 2 


w + 


^2,1 

#2,2 


+ 


£*21,0 

0 
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X = 


fn-p+m o 


V + 


Ip—m o 


I< 


rl 


y i 


= -F* 


Partitioning F along the lines of x: F = [F lt F 2 ], and with the usual definition of recovery error: 

E r (s) = T zu (s)M r (s), 


it works out that: 


with: 


M r (s) = F 2 (si — A r + K r C r ) 1 ( B r — h r D r ) 


A r — A 22 1 Ct — 


C'2,02 

^12 


£>, = 


F?21,0 

Si,i 


Minimizing M r is equivalent to finding a state feedback matrix of the transposed system to minimize 
the responses — the same form as for the full order estimator. 

A straightforward LQ/Ricatti solution yields the appropriate observer gain matrix which then 


e 

Recovery Error 

500 

0.9555 

250 

0.9007 

125 

0.9135 

62.5 

1.0012 


Aga 


well as plots of the yaw rate command 


rtga.n, looking at these values of the recovery error, _ 

response (which does not necessarily improve with decreasing e), the value of 250 was settled upon. 
Note that for less than this value, the recovery error does indeed start to rise from this point 
Consider that the recovery error presented in the table is no longer just the magnitude of the M(s) 
matrix. The eigenvalues of the closed loop system corresponding to the e = 250 design are shown 
in Table 7. Although the order of the system has been reduced by 10, the performance is somewhat 
dismal. More detail can be seen in the controller activity and recovery error singular value plots 
(see figures 19 and 20. One can see from these the large (and perhaps unnecessary) increase in 
activity for higher frequencies. This occurs even with the large values oft used. 

Thus, again we seek a way of frequency weighting the error in the transposed system used to 
generate the observer gains. A fair amount of trial and error led to using a corner of 0.2 rad/s on 
all criteria for this system. This is probably much lower than what intuition would lead one to. 
However, it is not very surprising in light of the experience with the frequency weighting on the 
full order observer. However, the recovery errors are slightly better. 


e 

Recovery Error 

0.5 

0.7437 

0.25 

0.6915 

0.125 

0.6987 

0.0625 

0.7421 


Still, we have relatively large values of t compared to a weighted full order system, though, in line 
with intuition, they are somewhat smaller than before. The recovery error for the e = 0.25 was best, 
though the improvement does not get one near the residual errors for a shaped full state system. 
In fact, it does not do better than an unweighted full order system. The corresponding eigenvalues 
for this system are shown in Table 8. Singular value plots of the controller and the recovery error 
are shown in figures 21 and 22. 

Overall there seems to be a fairly high price to going to this reduced order observer structure. 
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Max Singular Value 


Poles 


Damping 

Freq (rad/sec) | 

-7.83040 

± 

51.6865i 

0.150 

52.2762 

.-9.38043 

± 

51.2638i 

0.180 

52.1149 

-12.4905 

± 

41.0027i 

0.291 

42.8630 

.-6.11002 

± 

37.25741 

0.162 

37.7550 

-12.7108 

± 

25.28261 

0.449 

28.2979 

.-8.72514 

± 

24.60111 

0.334 

26.1025 

.-24.5081 

± 

2.91 104i 

0.993 

24.6804 

-24.5195 

± 

2.279731 

0.996 

24.6252 

-3.71934 

± 

19.6099i 

0.186 

19.9595 

.-4.20761 

± 

19.5369i 

0.211 

19.9849 

-19.3829 



1.000 

19.3829 

-17.5160 



1.000 

17.5160 

-6.90951 

± 

4.917291 

0.815 

8.48063 

.-5.98385 

± 

5 .242031 

0.752 

7.95521 

.-6.51433 



1.000 

6.51433 

.-3.75254 

± 

4.67419i 

0.626 

5.99413 

-5.79405 

± 

0.39518i 

0.998 

5.80751 

.-4.65048 



1.000 

4.65048 

.-1.87052 

± 

1.37995i 

0.805 

2.32446 

.-1.70573 

± 

0.10905i 

0.998 

1.70922 

.-1.48464 

± 

0.600511 

0.927 

1.60149 

.-4.26079e-3 

± 

6.65864e-3i 

0.539 

7.90518e-3 


Table 7: Reduced Order Observer Recovery, 2 Integral Control 



Frequency (rad/s) 


Figure 19: Recovery Error, Reduced Order Observer 
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Max Singular Value 



Poles 

Damping 

Freq (rad /sec) 

-8.80979 ± 52.3960i 

0.166 

53.1315 

-9.38043 ± 51.2638i 

0.180 

52.1149 

-16.1812 ± 35.1029i 

0.419 

38.6528 

-6.11002 ± 37 .257 4i 

0.162 

37.7550 

-14.0871 ± 24.0717i 

0.505 

27.8907 

-8.72514 ± 24.601H 

0.334 

26.1025 

-24.7999 ± 0.55088i 

1.000 

24.8060 

-24.5081 ± 2.91 104i 

0.993 

24.6804 

-20.5883 ± 1.420*241 

0.998 

20.6373 

-2.25247 ± 20.27891 

0.110 

20.4037 

-4.20761 ± 19.53691 

0.211 

19.9849 

-19.3829 

1.000 

19.3829 

-5.98385 ± 5. 2420 3 i 

0.752 

7.95521 

-6.54135 ± 0.58733i 

0.996 

6.56766 

-6.51433 

1.000 

6.51433 

-3.75254 ± 4.67419i 

0.626 

5.99413 

-4.72891 

1.000 

4.72891 

-4.65048 

1.000 

4.65048 

-1.87052 ± 1.379951 

0.805 

2.32446 

-1.70573 ± 0.10905i 

0.998 

1.70922 

-1.48464 i 0.6005H 

0.927 

1.60149 

-4.26079e-3 ± 6.65864e-3i 

0.539 

7.90518e-3 


Table 8: Reduced Order Observer Recovery, 2 Integral Control, Weighted 




Results for Various Weighting Factors 



Figure 21: Recovery Error, Shaped Reduced Order Observer 


Results for Various Weighting Factors 
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ieure 22: Controller Activity, Shaped Reduced Order Obser 




5.3 Numerical LTR I: 8th Order Design 


Aside from the structured observers, one can formulate a compensator with dynamics whose states 
do not necessarily have any analogues with observer states. This direct optimization run through 
the program SANDY [1], will optimize a controller design given a general structure tor it. In t is 
case, we specify that it. be 8th order in a series of 4 uncoupled 2nd order systems. The initial values 
for the direct feedthrough matrix on this controller correspond to those corresponding columns o 
the full state feedback matrix. The general control structure is as follows: 
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Of course the values marked with a 0 remain 0. The rest are allowed to optimize including the 
state feedback gains noted in the D c matrix. This still amounts to a hefty number of parameters: 
104 of them. It turns out that the feedthrough gains associated with states that turn out to be 
directly measured (i.e., pitch, roll, etc.) have a significant impact on the recovery. Optimization 
of these gains from their state feedback values provides for a significant amount of the recovery. 
Thus, an analogy from the Ricatti-based solutions does not hold here we are indeed changing the 
supposed state feedback gain matrix K. The recovery error is 0.16374, which is the best recovery 
of all systems. This is remarkable, considering the simple dynamics structure. Later on one will 
see the limits to optimizing the direct feedthrough gains associated with just the observed terms. 
This optimization perhaps could be improved upon as some analysis indicates that some further 
improvement could be made, however some current algorithmic difficulties have prevented it. The 
full state feedforward gains were used, as attempts to structure around their use did not seem to 
work (see figure 23). It would have seemed that recovery error could have been improved by the 
additional degrees of freedom, but the results seemed to indicate otherwise. This optimization 
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Zero 

Damping 

Freq (Hz) 

-9.37520 

± 

49.70025 

0.185 

50.5767 

-47.5601 



1.000 

47.5601 

-6.58311 

± 

37.7641i 

0.172 

38.3336 

-10.4636 

± 

27 20 17i 

0.359 

29.1448 

-20.6436 

± 

12.9692i 

0.847 

24.3795 

-23.4432 

± 

2.81882i 

0.993 

23.6121 

-3.74904 

± 

19.5247i 

0.189 

19.8814 

-11.2690 

± 

11.87101 

0.688 

16.3680 

-15.9453 



1.000 

15.9453 

-4.30385 

± 

6.15915i 

0.573 

7.51387 

-5.44168 

± 

4.74486i 

0.754 

7.21980 

-6.95416 

± 

1.78549i 

0.969 

7.17971 

-5.61014 



1.000 

5.61014 

-2.05374 

± 

1.28 164i 

0.848 

2.42084 

-1.63643 

± 

0.73012i 

0.913 

1.79192 

-1.50904 

± 

0 . 12381 i 

0.997 

1.51411 

-4.25281e-003 

± 

6.68477e-003i 

0.537 

7.92292e-003 


Table 9: Poles of 8 State LTR Optimized Controller and Plant 



Figure 23: LTR for Output Feedback Controller C(s), No Kjf (not used) 
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proceeded with filtered disturbance inputs. The pitch and roll command disturbance inputs were 
integrated, then passed through a filter with a time constant of 8 rad/sec. The vertical rate and 
yaw rate command inputs were passed through a filter with a time constant of 4 rad/sec. These 
filters were also applied in the subsequent SANDY based LTR optimizations. 

There is no control penalty R on an LTR-based direct optimization. Since it stands to reason 
that one is trying to recover a state feedback design, one needs to be able to freely vary the necessary 
gains. Ricatti-based LTR theory indicates a need for some gains to go to infinity, Fortunately 
most optimization schemes will detect an extreme lack of sensitivity of the cost function value 
to a parameter change: the parameter will probably never move very far before the optimization 
completes and terminates. Such has been the observed behavior here. 

The disturbance response penalty affected the mean-squared values. Note that here, as opposed 
to the closed form Ricatti solutions, one could minimize the actual error rather than a matrix 
analogous to some error multiplier matrix M(s). The value for the response penalty matiix Q, foi 
all direct LTR optimizations was: 

Q = diag{0. 1,0.05, 2, 1,1, 2, 1,0.2, 1,0.2} 

The entries corresponded to the outputs { p, q, r, u, v, w, <j), 0, f 0 , and f <p }. 

5.4 Numerical LTR II: 4th Order Design 

How much worse do the results get if the dynamics in the controller is lowered to 4t,h order? In a 
word, not much worse-a recovery error of 0.17105. With the overall structure otherwise the same, 
an interaction term between the pair of second order systems was allowed. This was intended to 


Zero 

Damping 

Freq (Hz) 

-464.184 

1.000 

464.184 

-8.95655 ± 51.1 120i 

0.173 

51.8908 

-6.90773 ± 37 .09391 

0.183 

37.7316 

-7.07229 ± 25.5506i 

0.267 

26.5114 

-24.1159 ± 2.94300i 

0.993 

24.2948 

-4.51612 ± I9.8662i 

0.222 

20.3731 

-20.2468 

1.000 

20.2468 

-5.75452 ± 5.30307 i 

0.735 

7.82541 

-2.82475 ± 5-709971 

0.443 

6.37048 

-4.12753 ± 5.42571i 

0.605 

6.81724 

-6.75926 

1.000 

6.75926 

-2.12814 ± 2.73388i 

0.614 

3.46455 

-1.38437 ± 1 . 1 3227 i 

0.774 

1.78844 

-1.72118 

1.000 

1.72118 

-1.57803 ± 0.61 632i 

0.931 

1.69411 

-1.12528 

1.000 

1.12528 

-4.2532e-003 ± 6.68453e-003i 

0.537 

7.92290e-003 


Table 10: Poles of 4 State LTR Optimized Controller and Plant 


make the dynamics a little richer to see if this would allow results nearly as good as those for the 
8th order system. Because the size of this controller was minimal it was felt that some extra effort 
needed to be made to see if this controller would achieve a level as good as any other system. These 
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interaction terms-from the output of one second order system to the input of the other made the 
design of the controller matrix look like: 

'0 1 0 O' 

_ A21 A22 ^23 0 

Ac - 0 0 0 1 

A4 1 0 >443 .444 

Since the cost (needing to optimize only 2 more variables) was low, this was done. 


5.5 Numerical LTR III: Oth Order Design 

Finally, one asks, is there a need for dynamics at all? What about feeding the output directly back? 
This time the direct feedthough matrix is all there is, so one will begin with using the state feedback 
gains and try to optimize it further. Analyzing results with just state feedback gains applied to 
the measured output, we see what one can expect in the least. The recovery error is 6.1807, with 
the accompanying suite of system eigenvalues as seen in table 11. While the performance is not 
horrible, it does not keep the response within any sort of acceptable tolerance. The eigenvalues of 
the system (Table 11) still indicate a well stabilized system. If one optimizes this D matrix via 


Zero 

Damping 

Freq (Hz) 

-9.49008 

± 

5 1 .5 146i 

0.181 

52.3815 

-5.48053 

± 

33.7356i 

0.160 

34.1778 

-7.55594 

± 

24. 5267 i 

0.294 

25.6642 

-24.3200 

± 

2.761441 

0.994 

24.4763 

-4.68312 

± 

20.9348i 

0.218 

21.4522 

-19.7090 



1.000 

19.7090 

-10.4723 



1.000 

10.4723 

-5.46663 

± 

5. 15352i 

0.728 

7.51284 

-3.40089 

± 

5 468781 

0.528 

6.44000 

-2.20360 

± 

2.526 14i 

0.657 

3.35220 

-1.80756 



1.000 

1.80756 

-1.47456 

± 

0.83193i 

0.871 

1.69305 

-1.34444 

± 

0.50772i 

0.936 

1.43711 

-4.10885-003 

± 

6.72672e-003i 

0.521 

1.25451e-003 


Table 11: Poles of UH-60 Stabilized With Selected State Feedback Gains 

SANDY, the recovery error improves dramatically to 1.2489. While this is still not as good as the 
optimized LTR designs with dynamics, it is not so far off of a Luenberger LTR design. Technically 
there are dynamics even with this design in the form of the integrators. 

Looking at the plot of recovery error maximum singular values (figure 24), one will see that the 
error at lower frequency is actually better for the Oth order controller. It is believed that this is 
a function of how well the optimization did. Because there were many fewer parameters to move, 
the relative strength of the optimization was probably better. Hopefully future work will improve 
technique to allow opitmization to the same degree on the 8th and 4th order designs. Because of 
the lack of internal dynamics, the error was worse than the 8th and 4th order designs for frequencies 
in the 1-10 rad /sec range. 
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Max Singular Value 


Zero 

Damping 

Freq (Hz) j 

-9.78679 ± 5 1 .42 lOi 

0.187 

8.33081 

-2.09033 ± 35.181 li 

0.059 

5.60912 

-7.35643 ± 24.4792i 

0.288 

4.06810 

-24.2160 ± 3 .31 6 15i 

0.991 

3.89007 

-5.59656 ± 19.3922i 

0.277 

3.21233 

-20.0660 

1.000 

3.19360 

-3.54972 ± 6.16242i 

0.499 

1.13186 

-5.86258 ± 5.180501 

0.749 

1.24515 

-6.71038 

1.000 

1.06799 

-2.87716 ± 4.32472i 

0.554 

0.82671 

-1.15676 ± 0.77973i 

0.829 

0.22202 

-1.10283 ± 0.67672i 

0.852 

0.20593 

-1.70616 

1.000 

0.27154 

-4.10885-003 ± 6.72672e-003i 

0.521 

1.25451e-003 


Table 12: Poles of 0 State LTR Optimized Controller and Plant 


Results for Various Order Controllers 
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Figure 24: Recovery Error, 8th, 4th and Oth Order 
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Figure 25: Controller Activity, 8th, 4th and Oth Order 


6 Direct Optimization for a 2 Integral State System Controller 


What is the difference between LTR and a direct optimization? Formwise, the compensator struc 
ture is identical to that of the 4th order LTR design: 
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The disturbance inputs, while filtered as before, will influence the dynamics in a significantly 
different way. The feedforward matrix normally used is based on the full state feedback design— 
these are the gains used to translate a vertical rate command into commands on the actuators. With 
this design, the yaw rate and vertical rate disturbances are summed into the signals going from the 
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sensors to the controller itself-no "recovery” to a set of full state feedback gains is implied (see 
figure 26). The pitch and roll commands, as before, go into the integrators (no direct feedthrough 



Figure 26: Direct Optimization for Output Feedback Controller C(s) 

gains around these integrators). Because there is no direct comparison with the original state 
feedback design feedforward, one cannot look at residuals between the two systems in a meaningful 
way. The overall recovery process appears somewhat sensitive to the feedforward gains. One can, 
however, look at the basic responses of this closed loop system and judge how well the optimization 
went. In fact, because there exists a 4th order LTR, design with the same structure a comparison 
to it will be useful. 

The largest source of difficulty was in choosing the relative weighting of control versus criterion. 
Classically, one worked for the fastest adequate response time versus having unreasonable gains or 
a large overshoot. The weighting matrices had values as shown in tables 13 and 14. Considering 


Q for State 
Feedback 

Q for Optimized 
Output Feedback 

Factor 

Weighted 

r i 

10 ~ " 

Blend of: 0.25«/ + 9 + f 0 dt 

0.1 

6 

Blend of: 0.2p + (fr + / <fi dt 

1.6 

16 

Yaw Rate r 

2 

30 

Vertical Rate w 


Table 13: Diagonal Values of Q and Corresponding Factors 


how this design had some high frequency components and overshoot, the difference in weights from 
the LQ to the direct design is significant. The resulting eigenvalues are shown in lable 15. From 
previous experience, one can see that the eigenvalues do not say much by themselves. Looking at 
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R for State 
Feedback 

R for Optimized 
Output Feedback 

Factor 

Weighted 

10 

0.4 

Collective 6q 

1 

1 

Sine 8 S 

1 

1 

Cosine 8 C 

1 

1 

Tail Rotor 


Table 14: Diagonal Values of R and Corresponding Factors 


Zero 

Damping 

Freq (Hz) 

-470.046 



1.000 

470.046 

-115.173 



1.000 

115.173 

-9.39606 

± 

5 1 .4281 i 

0.180 

52.2794 

-2.85548 

± 

35.6182i 

0.080 

35.7325 

-30.9068 



1.000 

30.9068 

-8.59788 

± 

26.2384i 

0.311 

27.6112 

-24.2882 

± 

3.56569i 

0.989 

24.5486 

-5.44497 

± 

19.9205i 

0.264 

20.6512 

-17.0312 



1.000 

17.0312 

-1.23462 

± 

8.68975i 

0.141 

8.77702 

-3.91763 

± 

6 . 1 27 17i 

0.539 

7.27255 

-2.88111 

± 

3.14284i 

0.676 

4.26359 

-3.07029 

± 

1.37917i 

0.912 

3.36583 

-3.74960 



1.000 

3.74960 

-1.34781 

± 

0.58055i 

0.918 

1.46753 

-1.46438 



1.000 

1.46438 

-0.33784 



1.000 

0.33784 

-3.96107e-003 

± 

6.69608e-003i 

0.509 

7.77995e-003 


Table 15: Poles of 4 State Direct Optimized Controller and Plant 
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Figure 27: Controller Activity, non-LTR Optimization 

the controller activity (see Figure 27), one can see that it starts out with a fairly large gain at low 
frequencies and increases a modest relative amount as one goes up in frequency. Considering the 
command response plots (figures 28 to 31, this extra controller activity is apparent in the faster 
response times (and overshoots) for pitch and roll along with the appearance of high frequency 
components in the responses. Other than this, there seem to be other differences with the LTR 
design, in that the direct optimization performed better in keeping the yaw rate and up rate near 
zero in the pitch and roll commands, as well as mitigating actuator use in the yaw rate and up rate 
commands. Note that while the rise time of the direct optimization response for the pitch and roll 
commands was better, the delay time was the same. 

The Bode plots do not add much to the insights already learned aside from reinforcing the 
fact that the rolloff is lower on the direct design. Because the phase also moves more slowly the 
difference in robustness may only be slight. 


7 Conclusions 

Though one could directly optimize a controller, the LTR procedure is in itself useful combined with 
the numerical optimization as the process derives a best case controller (versus full-state feedback) 
as well as a feel for a particular plant model behavior. Though Q and R will almost always be dealt 
with as diagonal matrices, developing the weightings for a proper controller design is one of the 
more difficult parts of the design — this development is sped up considerably given the availability 
of a full state design with which to explore the various effects. 

Recovering to a design (though it be to an unrealizable full state one in general) is well defined 
within a numerical LTR. Though the weighting matrix on the controls is zero, the gains do not 
tend to blow up (an optimization process would be hard pressed to do this). Having a full state 
design to recover to (as opposed to a direct optimization to the best possible controller) represents 
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Roll rate u r: Yaw rate Tail Rotor Roll Collective 



Responses to Pitch Command 



‘ 0 5 10 0 



Responses to Pitch Command 



Responses to Pitch Command 

Figure 28: LTR and Direct Comparison 








Roll rate u r: Yaw rate Tail Rotor Roll Collective 



Responses to Roll Command 



Responses to Roll Command 



Responses to Roll Command 

Figure 29: LTR and Direct Comparison 
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Roll rate u r: Yaw rate Tail Rotor Roll Collective 



Responses to Yaw Rate Command 



Responses to Yaw Rate Command 

Figure 30: LTR and Direct Comparison 
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Roll rate u r: Yaw rate Tail Rotor Roll Collective 



Responses to Up Command 



Responses to Up Command 



Responses to Up Command 

Figure 31: LTR and Direct Comparison 
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Figure 32: Pitch/Pitch Command Bode Plot, LTR and Direct 
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Figure 33: Roll/Roll Command Bode Plot, LTR and Direct 
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Frequency (Rad/Sec) 



Frequency (Rad/Sec) 

Figure 34: Yaw /Yaw Commanded Rate Bode Plot, LTR and Direct 
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Frequency (Rad/Sec) 


Figure 35: Up/Up Commanded Rate Bode Plot, LTR and Direct 
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METHOD 1 

ERROR 

Full State Observer (Unweighted) 

0.4359 

Full State Observer (Weighted) 

0.2746 

Reduced Order Observer (Unweighted) 

0.9007 

Reduced Order Observer (Weighted) 

0.6915 

8th Order Numerical LTR 

0.1637 

4th Order Numerical LTR 

0.1711 

0th Order Numerical LTR 

1.2489 


Table 16: Overall List of Recovery Error versus Method 


a continual design check— a sort of quality control. This is necessary as one may need to restart 
the optimizer to avoid some local minimum. In addition the LTR numerical solution may provide 
a good starting point for the direct optimization. 

Using the appropriate entries of the full state feedback gain matrix for the feedfoiward fiom the 
yaw rate command and up rate command to the actuators may be restrictive enough to affect the 
resulting controller. It appears that LTR does not work if these commands were inserted between 
the sensor output and the controller input as in the direct optimization (see Figure 26). Fortunately 
this is not an issue with the pitch and roll commands as the feed into the error integrators is the 
same in both LTR and the direct case. 

What would be the advantages and disadvantages of integral control on r and w? The process 
of recovery may actually be made easier; however, the full state feedback design would probably 
not have as good a response time as with a non-integral design. This sort of tradeoff may have to 
be made in making the decision as to what kind of mission one intends with this helicoptei. 

Other future exploration should include a test of the robustness of this controller to deviations 
about the hover flight condition. Ideally the gain scheduling from flight condition to flight condition 
should result in a fairly even performance across the envelope-especially at a midpoint between 
flight conditions. 

Further improvements to the numerical algorithms to combine robustness to defective matrices 
and the speed of the diagonal algorithms is also pending. Routine use of a robust form of optimiza- 
tion gradient calculation may improve the numerical results seen in this paper. The LTR process 
tries to recover a state feedback design with an output feedback design therefore it is inevitable 
that various eigenvalues, one set from the closed loop state feedback design, one set from the re- 
covery design, will overlap. Chances are this overlap will result in a defective set of eigenvalues. 
Currently this roadblock is overcome only by going to the much slower robust form. 
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